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PROCESSOR-SHARED  TIME-SHARING  MODELS  IN  HEAVY  TRAFFIC 

Donald  P.  Gaver 
Patricia  A.  Jacobs 


1 .   Introduction 

Processor  sharing  (PS)  is  a  mathematically  tractable  approxi- 
mation to  time  sharing,  a  procedure  followed  in  marr 
computer  systems.   In  effect.,  PS  assigns  to  each 
(i  =  1,2,...)  present  for  processing  1/ith  of  the  total  processing 
effort;  equivalently ,  a  single  30b  with  Markovian  service  rate 
u  completes  processing  in  (t,t+dt)  with  probability  (u/i)dt  +  o(dt) 
One  advantage  of  PS  is  that  short  jobs  are  not  tracked  behind 
long  jobs,  as  is  possible  in  a  FC-FS  discipline. 

Various  mathematical  results  have  been  obtained  aoout  certain 
processor  sharing  models.   An  early  example  was  the  oaper  of 
Coffman  et  al_.  (1970).   Recently  extensive  results  have  been  ob- 
tained for  Markovian  systems  by  D.  Mitra  (1981),  and  for  non- 
Markovian  single-server  Poisson  arrival  systems  by  T.  Ctt  (19  34) 
and  V.  Ramaswami  (1984). 

This  paper  is  a  continuation  of  work  reported  in  Gaver,  Jacobs, 
and  Latouche  (1984),  henceforth  GJL,  where  emphasis  was  placed  on 
proposing  and  evaluating  simple  approximations  to  the  distribution 
of  delay  experienced  by  a  particular  "tagged"  job  approaching  a 
time-shared  processor.   In  that  paper  it  was  shown  that  under  heavy 
traffic  conditions,  or  if  the  tagged  job  duration  became  large, 
then  the  distribution  of  taoqed  job  response  time  (also  called 
sojourn  time  by  others)  approaches  the  normal  or  Gaussian 
distribution . 


We  first  show  that  the  analysis  of  our  system  and  the  problem 
solution  can  naturally  and  conveniently  be  conducted  in  work  time 
rather  than  in  ordinary  clock  time.   Subsequently  we  focus  upon 
heavy-traffic  approximations  to  the  distribution  of  response  time, 
R(T) ,  given  the  actual  work  (computer  time)  requirement,  T. 
Specific  models  proposed  and  examined  are  (i)  for  a  system  of 
many  identical  terminals  that  independently  submit  jobs  or  pro- 
grams according  to  the  same  Markovian  process,  and  (ii)  for  a 
system  having  two  terminal  types,  each  of  which  submits  jobs  in 
a  manner  governed  by  its  own  Markov  process.   The  methodology 
extends  to  general  k  terminal  types  as  well,  and  to  other  models. 

The  approximation  solutions  are  evaluated  for  accuracy  by 
means  of  Monte  Carlo  simulations. 


2 .   The  Work-Time  Concept 

Imagine  that  a  tagged  job  requiring  T  units  of  processing 
approaches  the  computer.   Assume  that  it  arrives  when  the  system 
is  in  steady  state.   After  that  initial  moment  it  undergoes 
processing,  at  various  rates  governed  by  the  amount  of  its 
accompaniment,  until  T  units  of  service  or  work  are  accumulated, 
at  which  point  it  departs  after  a  random  delay  of  R(T) .   In 
following  the  tagged  job's  delay,  it  turns  out  to  be  convenient  to 
measure  time  in  terms  of  the  amount  of  actual  work  or  processing 
that  has  been  accomplished  on  the  tagged  job.   Thus  let  {X(w) , 
w    0}  denote  the  number  of  programs  or  jobs  undergoing  service 
at  a  moment  when  exactly  w  units  of  processing  have  been  accomplished 
on  the  tagged  job.   The  instantaneous  rate  of  accrual  of  clock  or 
response  time  at  work  time  w  is  clearly  X(w) :   if  X(w)  =  1  then 
the  tagged  job  is  alone  and  response  (clock)  time  and  work  time 
advance  at  the  same  rate,  while  if  X(w)  =  17  the  tagged  job  is 
accompanied  by  16  others  and  17  units  of  response  time  accrue  for 
every  single  work  time  unit.   It  follows  that 

T 
R(T)   =    /   X(w)  dw  .  (2.1) 

0 

For  the  models  to  be  considered,  the  process  X(w)  is  a  birth- 

and-death,  or  simple  Markov,  process  related  to  N(t),  the  number 

of  jobs  in  the  system  at  clock  time  t;  its  transition  rates, 
after  adjusting  for  tagged  job  entry,  are  seen  to  be 

Xfcdt   =   A(N-i)dt   =   A(N-i)idw   =   Awdw 


2.2) 


t  w 

y  .  dt   =   ydt   =   y  i  dw   =   y  .  dw 
l  l 


for  1    i    N  where  N  is  the  total  number  of  terminals  in  the 
system.   The  term  idw  is  required  to  allow  the  work-time  process 
to  advance  appropriately  in  clock  time.   Henceforth  we  drop 
the  superscripts,  allowing  the  context  to  imply  the  appropriate 
rate . 

The  approach  taken  here  invokes  the  work  time  concept 
described  above  to  facilitate  calculations,  first  for  a  single 
terminal  or  job  type  situation,  but  later  on  for  a  system  with 
more  diversified  traffic  patterns. 


3.   Heavy  Traffic  Analysis  of  a  Single  Terminal  Type  Markovian 
System 

Suppose  that  N  terminals  have  access  to  a  single  computer. 
Each  terminal  has  Markovian  demand  rate  X,  and  expected  service 
time  y   .   Service  times  are  assumed  to  be  independently  and 
exponentially  distributed.   The  discipline  at  the  computer  is 
PS.   It  is  clear  that  if  N(t)  is  the  number  of  terminals  that 
have  submitted  jobs  that  are  undergoing  service  at  the  computer 
at  (clock)  time  t,  then  {N(t) ,  t    0}  is  a  Markov  process  in 
continuous  time  that  is  identical  to  the  classical  single  repairman 
problem;  see  Feller  (1957),  p.  416.   This  is  so,  since  if 
N(t)  =  i  >  0,  then  each  individual  job  or  program  receives  (dt/i) 
units  of  processing  time  in  (t,t+dt) ,  and  hence  departs  with 
probability  y(dt/i)  +  o(dt) ,  but  the  probability  that  some  job 
departs  is  iy(dt/i)  +  o(dt)  =  ydt  +  o(dt) .   It  has  been  shown  by 
Iglehart  (19  65)  and  by  Burman  (19  79)  that  under  heavy  traffic  con- 
ditions (N  -*■   °°)  one  may  approximate  N(t)  by  a  suitable  Gaussian 
process,  namely  the  Ornstein-Uhlenbeck  process.   This  fact  alone 
enables  one  to  study  the  distribution  of  R(T) ,  and  to  deduce 
approximate  normality;  see  GJL  for  a  first  analysis. 

3 . 1   Diffusion  Approximation  in  Work  Time. 

Here  is  a  diffusion  process  approximation  for  X(w) .   On  the 
basis  of  intuition  write  down  the  stochastic  differential  equation 
for  X(w) ,  the  approximation  to  X(w) : 

dX(w)   =  X  [N-X(w)  JX(w) dw  -  yX(w)dw 

(3.1) 


+   />.  [N-X(w)  ]X(w)  +  yX(w)  dB(w) 


where  {B(w)  ,  w  >_  0 }  is  standard  Brownian  motion.   The  first 
right-hand-side  term  represents  infinitesimal  drift  of  X(w), 
while  the  second  is  the  diffusion  or  infinitesimal  variance  term, 
the  form  of  which  is  obtained  from  the  observation  that  arrival 
and  departure  processes  compete  like  independent  Poisson  processes 
in  short  time  periods.   Now  suppose  that,  as  N  -*-  «> , 

X(w)   -   Nm(w)  +  /N  Z(.w)  (3.2) 

where  m(w)  is  a  deterministic  function  of  time,  and  { Z (w) , 

w    0}  is  a  stochastic  process,  the  properties  of  which  must  be 

discovered.   Substitute  (3.2)  into  (3.1)  to  obtain 

Ndm(w) +/NdZ (w)   =   A  [N-Nm(w) -/NZ (w)  ]  [Nm(w) +/NZ (w)  ] dw 

-    p  [Nm(w)+/NZ (w)  ]dw  (3.3) 


+    A  [N-Nm(w)  -/NZ  (w)  ]  [Nm(w)  WNZ  (w)  ]+y  [Nm(w)  +  /NZ  (w)  ]    dB  (w) 

Next  isolate  terms  of  order  N  and  /N;  the  result  is,  after 
stipulating  that  X1  =  AN,  a  constant  as  N  ■>  oo , 


0(N):   dg^w)   =   a1 [l-m(w) Jm(w)  -  ym(w)  ,  (3.4) 


0(/N) :   dZ(w)   =   {A ' [l-2m(w) ]-y}Z(w) dw 


+  /A  *  [l-m(w)  Jm(w) +ym(w)  dB (w)  ;       (3.5) 

the  stochastic  differential  equation  (3.5)  is  of  Ornstein- 
Uhlenbeck  (O-U)  form;  see  Arnold  (1974). 


Next  obtain  the  approximate  long-run  mean  as  the  solution 
of  (3.4)  with  dm/dw  =  0,  examining  only  the  heavy-traffic 
situation  in  which  A'  >  y: 


m(oo)   =   1  -  Ji_  ,     x'  >  y  (3.6) 


1  -  -±- 
AN 


If  the  above  solution  is  used  to  define  the  stochastic 
differential  equation  parameters  there  results 

dZ(w)   =   (-  A  '+y)  Z(w)dw  +  /2y  (1  -  (y/\  ')  )  dB(w)  ,       (3.7) 

which  suggests  that  {Z(w)}  can  be  considered  an  O-U  process  with 
constant  coefficients;  namely 

dZ(w)   =   -pZ(w)dw  +  adB(w)  ,  (3.8) 

the  solution  to  which  is 

Z(w)   =   Z(0)e"pw  +  a   f   e"pudB(u)  .  (3.9) 


2 
The  parameters  p  =  (A'-  y)  and  a      =    2y(l-y/A' 


3 . 2   Response  Time  Evaluation 
Let 


T  T 

R(T)   =    /   X(w)dw   -    /  [Nm(w)  +  /NZ(w)]dw         (3.10 
0  0 


approximate  the  response  time;  in  this  approximation  R(T)  is 
normally  distributed  (Gaussian).   First, 


T 
E[R(T)]   -   E[R(T)]   -   N   /   m(w)dw   =   N  ( 1  -  A)  T        (3.11) 

0  A 


Second, 


T 
Var[R(T)]   -   Var[R(T)]   =   var[/N   /   Z (w) dw]  (3.12) 

0 


2 


T  T  T 

=   N{E[  /  Z(w)dw   /  Z(u)du)  ]-(E[  /  Z  (w)  dw])  ^ }  ; 
0  0  0 


for  ease  of  writing  we  have  left  the  initial  condition  Z(0) 
implicit.   In  order  to  evaluate  the  above,  recall  that  the  tagged 
job  approaches  the  server  when  the  latter  is  in  equilibrium, 
i.e.,  at  t  =  °° .   It  may  be  shown  that  the  diffusion  approximation 
for  N(t) ,  the  number  undergoing  service  at  clock  time  t,  is 

N(t)   =   Na(t)  +  /N  Y(t)  (3.13) 

where  a(t)  is  a  deterministic  function  of  time  and  (Y(t)}  is  a 
particular  Ornstein-Uhlenbeck  process.   A  similar  analysis  to 
that  leading  to  (3.6)  and  (3.7)  yields 

a(«)   =   1  -  ^j    ,   NA  >  u    ,  (3.14) 


1  -  -ii- 

1    X' 


and 


E[Y(«)]   =   0,    Var[Y(<=°)]   =  ^L      =   ^H_  .  (3.15) 


see  GJL  for  a  derivation;  (3.13)  provides  the  initial  condition 

for  evaluating  moments  of  R(T)  ,  using  (3.10)  .   Identify  Z(0)  , 

the  initial  value  of  the  work  time  noisp  process  Z  (w)  ,  with  Y(°°)  . 

T 
According  to  (3.9) ,  this  implies  that  E[  /  Z(w)dw]  =0.   In  order 

0 
to  compute  the  Var[R(T)]  it  is  next  necessary  to  evaluate  the 

following  integral: 


T  T 

I(T)   =   E[  /   Z(w)dw   /   Z(u)du] 
0  0 


T  T 

E[2   /   Z(w)dw   /   E[Z(.u)|Z(w)J  du] 
0  w 


T  T 

2E[  /   Z(w)dw   /   Z(w)e"p(u"w)du]  (3.16 

0  w 


T 

=   2E[  f   Z(w)dwZ(w)-(l  -e"p(T"w)  )  ] 


T 

=   2E[  /   E[Z(w)2|  Z(0)  ]-(l  -e"p(T_w)  )dw]  . 
0  P 

Square  (3.9)  and  take  the  expectation  to  see  that,  condition- 
ally on  Z(0) , 


I(T)        =       2E[    /T(Z(0)2e-2pw+^-(l-e-2pW))^(l-e-p(T-w))dw 
0  2p  p 


2E[Z(0)2][     fTe-2pw(l-e-p(T-w))dw] 
P  0 


+    4      /T(l-e"2pw) (l-e-p(T-w))dw 
P       0 


Now    put    E[Z(0)2]    =    E[Y(°°)2]    =    u/X  '    =    a2/2p    to    see    that 


I(T)       =      4      /T[l-e-p(T"w)]dw      =      4tT-^d  -e"pT)] 
p2   o;  p2  p 

2 
q    T  ri  1   ,,  -pT.  , 

=    — [1-^(1  -e  )] 


Thus    it    follows    that 

2 

Var[R(T)]       -      Var[R(t)]       =      NT   2_[ 1  - -i(  1  -  e"pT)  ]    .  (3.14) 

A  pT 

P 

To  terms  of  order  T  this  agrees  with  (4.10)  of  GJL;  not  surprisingly 
the  additional  factor  in  (3.14)  can  actually  provide  numerical 
results  superior  to  those  of  GJL. 

The  form  of  the  heavy  traffic  approximation,  namely  the  limit- 
ing normal  form  with  parameters  (3.11)  and  (3.14) ,  can  be  more 
rigorously  validated  by  use  of  the  theory  of  convergence  of 
suitably  normalized  sequences  of  semigroups  of  transformations; 
see  Burman  (19  79) .   Details  appear  in  an  appendix  to  this  paper. 


L0 


4 .   Heavy  Traffic  Analysis  of  a  K-Terminal-Type  Processor 
Sharing  System. 

Consider  the  following  natural  extension  of  the  previous 

model.   The  processor  is  jointly  utilized  by  K  sets  of  terminals, 

each  generating  distinctive  job  types.   There  are  N.  terminals 

in  the  ith  set,  and  arrival  rate  and  service  rate  are  A.  and  u. 
—  1       i 

respectively.   Again  the  discipline  at  the  computer  is  PS.   Of 
course  this  is  not  the  same  as  a  situation  in  which  all  terminals 
are  the  same,  but  Type  j  jobs  occur  with  probability  p.  from  each 
terminal.   The  latter  model  can,  however,  be  studied  in  an 
analogous  heavy-traffic  manner,  as  can  other  interesting  models. 

4.1   A  Diffusion  Model  for  the  Work-Time  Process. 

Let  (X.(w),  i  =  1 , . . . , K}  represent  the  number  of  jobs  of  all 
types  present  at  the  computer  at  work  time  w.  The  present  model 
implies  that  {X. (w) }  is  a  multivariate  or  vector-state  birth  and 
death  Markov  process.  We  choose  to  study  a  diffusion  approxima- 
tion {X. (w) }  to  {X. (w) }  that  is  described  by  the  following  system 
of  s .d.e . : 


K  - 

dX.  (w)   =   X  .  (N.  -X.  (w)  )  (  I    X  (w)  )  dw  -  u  .  X.  (w)  dw 
i  ill     i,  _  -i  ^  ii 


A  .  (N.  -X.  (w)  )  I    X  (w)  +u-X.  (w)  dB.  (w)  (4.1! 

ill    ^  =-,  k       l  i       i 


i  =  1,2 ,K 


where  {B.(w)l  are  mutually  independent  standard  Brownian  motion 


1  1 


or  Wiener  processes.   The  work  time  process  is  a  transformation 
of  the  clock  time  process;  in  particular,  the  drift  of  the  ith 
component  of  the  clock  time  process  {N. (t)}  is  seen  to  be 

N. (t)dt 
Ai(Ni  -Ni(t)  )dt  -  Ui  -j- ,  (4.2; 

I       N  (t) 
k=l   K 


which  exhibits  the  processor-sharing  effect  in  the  term  multiply- 
ing i_i  .  .   Multiplication  by  the  total  in  service,  ^N.(t),  converts 
to  the  work  time  transition  rates,  in  analogy  with  (2.2)  . 
Now  once  again  approximate  by  writing 


Xi(w)   =   Nimi(w)  +  /Ni  Zi(w)  ,    i  =  1,2,. ..,k;       (4.3) 


m. (w)  and  {Z. (w) }  are  to  be  determined,  subject  to  the  normaliza- 
1         K  x 


tion  N  =     N,  ■*   °°  but  with 

kii^ 


N./N   ■*■  I.     ,         0  <  i .     <    1  ,  (4.4a) 

l         l        —   l  — 


and 


NA.   ■*   A  !  ,   0  <  A  !  <  °°,   A  .'  >  y  .  .  (4.4b) 

i       l-i-ii 


Conditions  (4.3)  and  (4.4)  are  referred  to  as  the  heavy  traffic 
normalization  (HTN) .   The  result  of  isolating  terms  according  to 
order  in  (4.1)  is : 


12 


dm.(w)  K 

0(N):   "" dw" —   =  \'i(l-mi(u))     I    ^kmk(w)  -uimi(w)  (4.5) 


k=l 


0(/N) :   dZ 


K  K 


(w)   =   A!{-(  \%  m,  (w)  )  Z  .  (w)+/T~(l-m.  (w)  )  Y  /I7z,(w)}dw 
1  X         k=lK  K      x  x  x  k=l   k  k 


-y.Z.(w)dw  (4.6) 


K 
A|(l-mi(w))  J  £Rmk(w) +yimi (w)  dBi ( 


for  i  =  1,2, ... ,K.   Thus  (m.  (w)  ;  w  >_  0;  i  =  1,2,  ...  ,K)  must 
be  found  by  solving  a  system  of  ordinary  first-order,  but  non- 
linear differential  equations,  while  (4.6)  shows  that 
{Z.(w);  w    0,  i  =  1,2, ...,K}  is  a  multivariate  Ornstein- 
Uhlenbeck  process. 

4.2   A  Diffusion  Model  for  the  Clock  Time  Process. 

In  order  to  provide  the  initial  conditions  encountered  by 
the  tagged  job,  it  is  necessary  to  study  the  clock-time  process 
N.(t);  see  (3.13).   The  corresponding  approximation  has  s.d.e. 


N.  (t) 
dN.(t)        =      A.(N.-N.(t))dt-y.    -=-± dt 

1  111  1         Ja 

I     N     (t) 
k=l    K 

(4.7 


/  ~  N.(t) 

+    \/X  .  (N.-N.  (t)  )+y  .— dB.(t)        i    =    1,2,  ...,K 

ylll  l    i\  l 

I       N     (t) 
k=l      K 


1  5 


Now    invoke    the    HTN : 


N. (t)       =      N.a. (t)     +    /N.    Y. (t)     , 
1  11  11 


(4.8) 


K 
and   again   N   =  N,     ->■   °° ,    with 

k=l      K 


N./N      ■*      i  ■     ,         0    <    £  .    <    1    , 
l  l  —      l    — 


4.9) 


but 


yi/N 


y !    ,        0    <   y .'    <   °° ,      X  .    >    y  !  . 

i  —     i  li 


4.10 


The  result  of  isolating  terms  is 


0(N] 


dai(t: 
dt- 


X.(l.a.(t))  -  y!  ^ 


aj_(t) 


I  ^av(t) 


k=l 


k"k 


(4.11) 


0(/N) 


dYi(t; 


-XiYi(t)dt-u-{- 


Yi(t) 


I    ^vav(t 


k=l 


k  k 


I    ^vYv(fc 


+  /FTa.  (t) 


k=l 


(  I  4.  a.  (t)) 


}dt 


k=l 


k  k 


/                 ai(t) 
+  Wx1(l-ai(t))+u[  -^ — dBi(t: 


I  Vv(t) 


k=l 


k~k 


4  .12: 


These  equations  closely  resemble  those  describing  the  work  time 
approximation;  again  the  semigroup  approach  is  applicable. 

If  a  long-run  solution  to  the  0(N)  term  exists  in  work  time, 
and  consequently  dnK/dw  -*  0  as  w  -»■  °° ,  the  result  is  the  system 
of  equations  for  m.  (<»)  =  m.  : 

K 
pi(l-mi)  I    2,kmk  -  m..   =   0  ,  (4.13) 

K=  x 

where  p.  =  XVy.  =  NA  .  /y  .  .   Now  these  same  equations  are  satis- 
fied by  a  presumed  long-run  solution  in  clock  time,  i.e.,  if 
da./dt  ■*    0  in  (4.11);  for  a.(»)  =  a.: 


K 

I 
k=l 


p.(l  -ax)  -  a±/  I      ikak      =      0  .  (4.14 


Consequently  the  long-run  solutions  in  work  and  clock  time  agree 
at  the  0(N)  term  level;  this  means  that  the  long-run  mean  number 
present  in  both  clock  and  work  time  agree: 


E[Ni(t)]   -   Nai(-)   =   Nm^co)   -   EtX^t)]  .     (4.15) 


Next  substitute  these  long-run  results  in  the  s.d.e.  to  see 

that  as  t,w  ■+   °°,  Y.  and  Z.  are  essentially  the  same  process. 

11 
K 

Put    S    =      [    lvav    to    simplify   writing.       Then 

k=l    K    K 

Y.(t)  _   J/Vk(t) 

dY.  (t)       =      -X.Y.  (t)dt  -  y  '    -=^ dt    +    A  .  (1-a.  )  /?.  .    = dt 

1  ll  lS  ill  s 


+    /2A .  (1-a.  )     dB .     ,  (4.16) 


L5 


or 


dY,  (t)  =A.{-Y.  (e)  S  +/r"(l-a.)    /ITy,  ( t)  }^  -  p  !  Y  .  ( t)  ~ 

1  11  1       l,,i^K        O       11        O 

k=l 


1 


+  /2A.(l-a.)S  —  dB.  (t)  ,  (4.17) 

1     X    /S    X 


and  a  direct  comparison  with  the  corresponding  equation  in  work 
time,  (  4.6  ) ,    shows  that  the  long-run  behaviors  of  the  two 
processes  { Z . (w) }  and  {Y . ( t) }  are  identical  except  for  a  constant 
time-scale  change:   for  large  w  and  t, 


(Zi(w)}     and     {y..^)}  (4.18) 


have  the  same  probability  law;  i.e., 
finite-dimensional  distributions  and 
limiting  distribution. 


4.3  Response  Time 

We  discuss  the  response  time  under  these  conditions:   a 
tagged  job  approaches  the  processor  when  the  latter  has  been 
operating  for  some  time,  so  the  long-run  clock  time  distribution 
prevails;  after  arrival,  the  job  remains  present  until  the  total 
work  time  accumulated  on  the  job  is  T,  the  requested  service 
time,  giving 


lb 


T   K  K     T 

R(T)       /    I     [Xi(w)dw]  ~         I         j      X.(w)dw  (4.19) 


i=l  i=l  0 


K        T  K        T 

=  N  I       I.       j      m.(w)dw  +  /N  I    /l~.       j      Z.(w)dw 
i=l   x  0     x  i=l  x   0     x 


where  it  is  understood  that  the  initial  condition  for  the 
Z. (w)  integrand  in  (4.19)  is  given  by  the  approximate  stationary 
distribution  from  the  clock  time  process.   In  view  of  (4.18) , 
this  is  equivalent  to  removing  the  initial  condition  by  the 
long-run  distribution  of  the  work  time  process  itself. 

Since  the  long-run  situation  is  being  discussed  it  is  first 
necessary  to  solve  the  steady-state  version  of  (4.5) 


K 

0   =   A'(l-m.)   T  l.m,       -   y-m.,   i  =  1 ,  2  ,  .  .  .  ,  K.  (4  .  20  ) 

1        1,-.K.K         11 

k=l 


Then  the  solution  provides  parameters  for  the  long-run  version 
of  (4.6);  here  written  in  matrix  form 

dZ(w)   =   A  Z_(w)dw  +  a  dB  (4.21) 

Now  to  find  the  variance  of  R(T) ,  append  the  row 


K 

dZ,  , (w)   =    >  /T~Z.(w)dw  (4.22) 

k+1         •  -i   i   i 
i  =  l 
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to  the  former  drift  matrix  A  of  (4.21) ,  and  consider  the  system 

dZ_*(w)   =   A*Z*(w)dw  +  a*dB*  (4.23) 

the  solution  to  which  can  be  formally  written  out  in  terms  of  the 
appropriate  fundamental  matrix,  and  computed  in  terms  of  eigen- 
values and  eigenvectors  of  the  matrix  A*.   See  e.g.  Arnold 
(1974),  Chapter  8  and  Coddington  and  Levinson  (1955)  for  details. 
A  convenient  way  of  formalizing  the  calculations  is  actually  by 
using  Laplace  transforms.   Unfortunately,  no  truly  simple 
formulas  result.   Finally,  the  covariance  matrix  C(w)  of  the 
components  of  Z*(w)  satisfies  the  matrix  differential  equation 


dC(w) 

— ^-  =  (A*)C(w)  +  C(w)  (A*)  '  +  (a)  (a)  ' 


where  '  denotes  transpose;  the  initial  conditions  are  provided 
by  the  long-run  distribution  in  clock  time,  or  in  view  of  (4.18)  , 
of  the  work  time  process  {Z_}  itself.   it  is  the  K+l—  diagonal 
element  of  C (w) ,  evaluated  at  w  =  T  and  multiplied  by  N  that 
provides  the  required  approximate  Var[R(T)]. 


5.   Simulation  Studies  of  the  Accuracy  of  the  Normal  Approxima- 
tions to  the  Distribution  of  Response  Time 

In  this  section  we  use  simulation  to  study  the  numerical 
accuracy  of  normal  approximations  to  the  distribution  of  the 
response  time.   Two  continuous  time  Markov  chain  models  were 
simulated.   In  one  there  is  a  single  terminal  type;  in  the 
second,  a  two-terminal  type  system  is  examined.   Two  normal 
approximations  were  evaluated:   one  results  from  a  central 
limit  theorem,  and  the  other  results  from  applying  the  previously 
derived  diffusion  approximation  to  the  Markov  processes. 

5.1   A  Single  Terminal  Type  Markovian  System. 

Let  X(w)  denote  the  number  of  other  jobs  undergoing  service  at 
a  moment  when  exactly  w  units  of  processing  has  been  accom- 
plished on  the  tagged  job  for  the  single  terminal  type 
Markovian  model  of  Section  3.   Since  {X(w);  w    0}  is  a  Markov 
process  and 

T 
R(T)   =    /   (X(w)  +1)  dw  , 
0 

it  follows  that  there  are  constants  m(c)  and  o(c)  such  that 

R(T)  -  m(c) T 
a(c)  /f 

converges  in  distribution  to  a  standard  normal  distribution 
as  T  -  »  (Cf.  Keilson  (1979),  p.  121).  Call  this  a  central 
limit  theorem  (CLT)  for  such  a  process.   In  this  case 


L9 


m 


(C)    =    1  +  I      7T  (j)  j 


where  it  is  the  stationary  distribution  of  {X(w);  w    0}  and 
a  formula  for  evaluating  a(c)  is  given  in  Keilson  (1979) .   The 

CLT  normal  approximation  states  that  R(T)  has  a  normal  dis- 

2 

tribution  with  mean  m(c)T  and  variance  a {c)    T.   The  CLT  mean 

is  the  true  mean  for  R(T)  under  steady  state  (cf .  GJL) . 
The  CLT  normal  approximation  should  be  increasingly  accurate 
as  T  becomes  large,  despite  values  of  other  system  parameters, 
including  the  number  of  terminals . 

The  derivation  of  the  heavy  traffic  (or  diffusion)  approxi- 
mation is  detailed  in  Section  3.   In  summary,  the  HT  approxi- 
mation is  that  R(T)  has  a  normal  distribution  with  mean 


N ( 1  -  — ) T       JL  <  i 
m     XN;   '     AN    l 


and  variance 


2 
NT  2_[1  -  -^(1  -  e"pT)  ]  , 

P 


where 


a2      =   2ud  -  &> 


and 


p   =   NX  +  u 


The  mean  and  variance  of  the  HT  approximation  are  easier  to 
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compute  than  those  for  the  CLT .   It  is  anticipated  that  the  HT 
approximation  should  be  increasingly  accurate  as  N  becomes  large 
when  heavy  traffic  conditions  prevail,  i.e.   y^  <  1 .   It  is 
inapplicable  under  other  circumstances.   We  have  conducted  simu- 
lations to  assess  these  anticipations.   All  simulations  were 
carried  out  on  an  IBM  3033  computer  at  the  Naval  Postgraduate 
School  using  the  LLRANDOMII  random  number  generating  package 
(see  Lewis  and  Uribe  (1981) )  . 

Conditional  response  times  given  the  number  of  jobs  being 
processed  at  the  time  of  arrival  of  the  tagged  job  were  simu- 
lated; the  tagged  job  required  T  time  units  of  processing.   For 
each  initial  condition,  500  replications  were  done.   Sample 
moments  and  relative  frequencies  were  computed  for  each  initial 
condition  giving  conditional  response  time  sample  moments,  and 
selected  response  time  relative  frequencies,  i.e.,  estimated 
probabilities  of  response  times  in  selected  ranges.   Uncondi- 
tional sample  moments  and  relative  frequencies  were  then  computed 
by  multiplying  each  conditional  moment  or  relative  frequency 
by  the  appropriate  stationary  probability  of  there  being  j  jobs 
present  at  the  time  of  arrival  of  the  tagged  job  and  then 
summing  over  all  possible  j .   The  stationary  probability  is 
of  the  form  k A tt  ( j )  where  k  is  chosen  so  that  the  probabilities 
sum  to  1  (cf.  Kelly  (1979)).   A  detailed  description  of  the 
simulation  program  can  be  found  in  Pornsuriya  [1984]  . 

Simulated  and  approximating  means  and  standard  deviations 
for  various  values  of  N,  A,  and  y  appear  in  Tables  1-3.  Some 
discussion  of  specific  cases  now  follows.   When  N  =  10 ,  A  =  25 


i\ 


and  u  =  100,  it  follows  from  (3.6)  that  approximately 
(10)  (1  -  yen-)  =  6  jobs  are  being  processed  along  with  the 
tagged  job;  thus  the  traffic  is  moderate  in  this  case.   When 
N  =  10,  A  =  15,  \i    =  100,  then  on  the  average  3.3  are  being 
processed  with  the  tagged  job;  a  rather  light  traffic  case.   When 
N  =  25,  A  =  10,  y  =  100,  then  on  the  average  15  jobs  are  being 
processed;  again  a  moderate  traffic  case.   The  HT  mean  is 
lower  than  the  simulated  mean  for  N  =  10 .   For  N  =  25  it 
agrees  with  the  simulated  mean.   As  mentioned  before,  the  CLT  mean 
equals  the  true  mean.   The  CLT  standard  deviation  approaches 
the  simulation  value   as  T  becomes  large,  as  anticipated.   Also 
as  anticipated,  the  HT  standard  deviation  is  closer  to  the 
simulated  one  for  the  moderate  traffic  cases  than  for  the 
light  traffic  case.   In  order  to  assess  the  degree  of  normality 
of  the  distribution  of  R(T) ,  the  a-quantiles  for  each  approxi- 
mating normal  distribution  were  computed.   The  relative 
frequency  of  being  less  than  or  equal  to  each  a-quantile  was 
computed  using  the  simulated  data.   The  results  appear  in 
Tables  4-6.   For  N  =  25 ,  A  =  10 ,  y  =  100,  the  HT  approximation 
appears  to  describe  the  data  well  for  all  values  of  T.   For 
N  =  10,  the  HT  approximation  does  better  than  the  CLT  for  the 
moderate  traffic  case  of  A  =  25.   For  the  light  traffic  case, 
N  =  10,  A  =  15,  both  approximations  do  poorly  for  small  T  =  0.01, 
the  mean  work  request  time.   The  CLT  does  well  for 
large  T  =  0.10,  which  is  10  times  the  average  service  time. 

Note,  however,  that  all  simulations  have  been  carried  out 
for  modest  system  sizes,  N.   If  N  grows  to  say  50,  or  100,  the  HT 
approximations  can  be  expected  to  improve  correspondingly; 
they  are  often  not  bad  even  at  the  level  of  N  =  25. 
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TABLE  1 

Simulated  Mean  and  Standard  Deviation  for  R(T] 
and  Their  Approximating  Values 


N  =  10,  A  =  15,  \i    =    100 


TIME  T 

Mean 

Std.  Dev. 

0.01 

Simulation 

.0404 

( .0002) * 

.0172 

CLT 

.0404 

.0323 

HT 

.0333 

.0238 

0.025 

Simulation 

.1005 
(  .0005) 

.0375 

CLT 

.1010 

.0510 

HT 

.0833 

.0535 

0.05 

Simulation 

.2016 
(  .0010) 

.0622 

CLT 

.2019 

.0722 

HT 

.1667 

.0919 

0.10 

Simulation 

.4029 
(  .0015) 

.0940 

CLT 

.4039 

.1020 

HT 

.3333 

.1462 

The  standard  error  of  the  estimate  of  the  mean 


.!  \ 


TABLE  2 

Simulated  Mean  and  Standard  Deviation  for  R(T) 
and  Their  Approximating  Values 


N  =  10,  A  =  25,  u  =  100 


TIME  T 

0.01 

Simulation 

CLT 

HT 

0.025 

Simulation 

CLT 

HT 

0  .05 

Simulation 

CLT 

HT 

0  .10 

Simulation 

CLT 

HT 

Mean 

Std.  Dev. 

.0606 

.0503 

( .0002) * 

.0605 

.0245 

.0600 

.0160 

.1507 

.0315 

(  .0004) 

.1513 

.0288 

.1500 

.0314 

.3021 

.0497 

(  .0008) 

.3027 

.0548 

.3000 

.0481 

.6036 

.0748 

(  .0012) 

.6053 

.0776 

.6000 

.0706 

Standard  Error  for  the  Estimate  of  the  mean 
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TABLE  3 

Simulated  Mean  and  Standard  Deviation  for  R(t; 
and  Their  Approximating  Values 


N  =  25,  ,\  =  10,  y  =  100 


TIME  T 

0.01 

Simulation 

CLT 

HT 

0.025 

Simulation 

CLT 

HT 

0.05 

Simulation 

CLT 

HT 

0.0625 

Simulation 

CLT 

HT 

Mean 

Std.  Dev. 

.1498 

.0812 

.0002) * 

.1500 

.0390 

.1500 

.0254 

.3759 

.0503 

.0006) 

.3750 

.0617 

.3750 

.0497 

.7506 

.0804 

.0010) 

.7500 

.0872 

.7500 

.0760 

.9381 

.0909 

.0011) 

.9375 

.0975 

.9375 

.0863 

The  standard  error  of  the  estimate  of  the  mean 
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TABLE  4 

Simulated  Probability  (Relative  Frequency)  that  the 
Response  Time  is  Less  than  or  Equal  to  the  Approximating 

a-Quantile 


N  =  10,  A  =  15,  y  =  100 

TIME     a:      .10     .25    .50    .75    .90    .95    .99 
0.01 


CLT 

0.00 

.12 

.50 

.89 

1.0 

1.0 

1.0 

HT 

0.00 

.10 

.36 

.68 

.91 

.97 

1.0 

0.025    CLT     0.04    .21    .49    .81     .97    .99   1.0 
HT      0.00    .09    .33    .67     .91    .98   1.0 


0.05 


CLT 

0.08 

.23 

.49 

.77 

.94 

.98 

1.0 

HT 

0.0 

.07 

.30 

.65 

.92 

.97 

1.0 

0  .10 


CLT 

0.10 

.24 

.48 

.76 

.92 

.97 

1.0 

HT 

0.00 

.04 

.24 

.60 

.89 

.97 

1.0 

2  b 


TABLE  5 

Simulated  Probability  (Relative  Frequency)  that  the 
Response  Time  is  Less  than  or  Equal  to  the  Approximating 

a-Quantile 


N  =  10,  A  =  25,  u  =  100 

TIME      a:     .10  .25    .50    .75    .90  .95  .99 

0.01      CLT    .04  .15    .45    .85    1.0  1.0  1.0 

HT     .11  .22    .43    .72     .91  .98  1.0 


0.025 


0.05 


0  .10 


CLT 

.08 

.19 

.45 

.80 

.98 

1.0 

1.0 

HT 

.10 

.22 

.44 

.73 

.92 

.98 

1.0 

CLT 

.09 

.22 

.45 

.77 

.95 

.99 

1.0 

HT 

.11 

.22 

.44 

.71 

.91 

.97 

1.0 

CLT 

.10 

.24 

.47 

.76 

.93 

.98 

1.0 

HT 

.11 

.23 

.45 

.70 

.89 

.95 

1.0 
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TABLE  6 

Simulated  Probability  (Relative  Frequency)  that  the 
Response  Time  is  Less  than  or  Equal  to  the  Approximating 

ct-Quantile 


N  =  25,  X    =  10,  u    =  100 
TIME      a:      .10    .25    .50    .75    .90    .95    .99 
0.01 


CLT 

.03 

.15 

.46 

.86 

.99 

1.0 

1.0 

HT 

.11 

.23 

.46 

.73 

.91 

.97 

1.0 

0.025 


CLT 

.07 

.19 

.46 

.80 

.96 

.99 

1.0 

HT 

.11 

.23 

.46 

.73 

.91 

.97 

1.0 

0.05 


CLT 

.09 

.21 

.46 

.78 

.94 

.98 

1.0 

HT 

.11 

.24 

.46 

.73 

.91 

.96 

1.0 

0.0  62  5    CLT 
HT 


0  9 

.22 

.47 

.77 

.93 

.98 

1.0 

11 

.24 

.47 

.73 

.90 

.96 

1.0 
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5.2   Simulation  Results  for  Markovian  Model  with  Two-Terminal 
Types . 

In  this  subsection  we  describe  the  results  of  a  simulation 

of  the  general  K-type  Markovian  model  of  Section  4,  in  the 

case  of  K  =  2  sets  of  terminals.   As  in  Section  4  let  X. (w) 

represent  the  number  of  other  jobs  of  type  i  being  processed  when 

the  tagged  job  has  acquired  exactly  w  units  of  processing. 

As  before  the  response  time  for  the  tagged  job  requiring  T 

units  of  work  is 


T 
R(T)   =    /   [X,  (w)  +  X  (w)  +  1]  dw  . 
0 


The  process  {  (X,  (w)  ,  X_  (w))  ;  w    0  }  is  Markovian.   Hence 
again  R(T)  satisfies  a  central  limit  theorem  as  T  ->-  °°. 
The  normal  approximation  for  the  distribution  of  R(T)  result- 
ing from  the  central  limit  theorem  will  again  be  referred  to  as 
CLT. 

The  mean  term  (m,  + m„ ) T  for  the  heavy  traffic  approximation 
was  computed  by  solving  the  system  of  equations  (4.20)  for 
m,  and  m_ .   The  variance  term  for  the  approximation  was  computed 
by  solving  the  system  of  stochastic  differential  equations  (4.23) 
as  detailed  in  Arnold  [Corollary  (8.2.4)].   The  fundamental 
matrix  (Arnold  [p.  129])  was  found  by  computing  Laplace 
transforms  of  the  system  of  defining  differential  equations 
and  then  inverting  the  solution.   The  approximating  variance 
term  was  found  by  computing  the  variance  of  the  solution  of 
the  s.d.e.   As  in  the  case  of  one-terminal  type,  it  is  a 
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linear  combination  of  exponentials  and.  constant  terms.   Its 
exact  form  is  uninf ormative  and  will  not  be  given  here. 

Conditional  response  times,  given  the  number  of  jobs  of 
each  type  being  processed  at  the  time  of  arrival  of  the  tagged 
job,  were  simulated.   The  tagged  job  was  always  taken  to  be 
a  Type  1  job.   For  each  initial  condition,  300  applications 
were  carried  out.   Sample  moments  and  probabilities  (relative 
frequencies)  were  computed  for  each  initial  condition  giving 
conditional  sample  moments  and  probabilities  (relative  fre- 
quencies) .   Unconditional  sample  moments  and  probabilities 
were  computed  in  a  similar  manner  to  that  in  the  one-terminal 
type  simulation;  see  Pornsuriya  [1984]. 

Values  of  the  simulated  means  and  standard  deviations  and 
their  approximating  values  for  R(T)  for  various  cases  in  which 
N,  =  5  and  N~  =  5  appear  in  Tables  7-10 .   Once  again  the  CLT 
mean  is  the  true  steady  state  mean  for  R(T) .   The  means  and 
standard  deviations  of  R(T)  for  each  T  differ  surprisingly 
little  for  the  four  cases.   This  suggests  that  perhaps  the 
two-type  terminal  model  can  often  be  satisfactorily  approxi- 
mated by  a  one-type  model  in  which  the  arrival  rate  and  service 
rates  are  the  average  arrival  and  service  rates  in  the  two-type 
model.   Values  of  the  simulated  means  and  standard  deviations 
and  their  approximating  values  for  the  approximate  one-type 
model  with  N  =  10  ,  A  =  25  and  \i    =  75  appear  in  Table  11.   The 
values  for  the  approximate  one-type  model  are  acceptably 
close  to  those  for  the  two-type  model.   Note  that  the  quality 


of  the  HT  approximation  is  generally  quite  good,  even  though 
the  system  sizes  N,  and  N2  can  hardly  be  called  "large." 
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TABLE  7 

Simulated  Means  and  Standard  Deviations  for  R(T) 
and  Their  Approximating  Values 


Nl  =  5'  N2  =  5/  Al  =  20 '    X2    =  30,  ul  =  50,  M2  =  10° 


TIME  T 

0.01 

Simulation 

CLT 

HT 

0.025 

Simulation 

CLT 

HT 

0.0375 

Simulation 

CLT 

HT 

0.050 

Simulation 

CLT 

HT 

0 .0625 

Simulation 

CLT 

HT 

Mean 

0  .0713 
(0.0001) * 

0.0707 

0.0710 

0 .1783 
(0.0003) 

0.1766 

0.1775 

0.2664 
(0.0005) 

0.2650 

0.2663 

0.3570 
(0.0005) 

0.3533 

0.3550 

0.4439 
(0.0006) 

0 .4416 

0.4438 


Std.  Dev 
0.0135 

0.0204 
0.0132 
0 .0263 

0.0322 
0.0253 
0.0353 

0.0395 
0.0325 
0 .0414 

0.0456 
0.0384 
0 .0478 

0 .0510 
0.0435 


Standard  error  of  the  estimate  of  the  mean 
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TABLE  8 

Simulated  Means  and  Standard  Deviations  for  R(T] 
and  Their  Approximating  Values 


Nl  =  5'  N2  =  5'  Al  =  30'  X2  =  2°'  yl  =  100'  y2  =  50 


TIME  T 


Mean 


Std.  Dev 


0.01 


Simulation 


CLT 
HT 
0.025      Simulation 

CLT 
HT 
0.0375     Simulation 

CLT 
HT 
0.050      Simulation 

CLT 
HT 
0.0625     Simulation 

CLT 

HT 


0.0710 
(0.0001) * 

0.0716 

0.0710 

0.1777 
(0.0003) 

0.1789 

0.1775 

0.2672 
(0.0004) 

0.2684 

0.2663 

0.3567 
(0.0005) 

0 .3578 

0.3550 

0  .4461 
(0.0006) 

0 .4473 

0  .4438 


0.0137 

0.0204 
0.0132 
0.0268 

0 .0322 
0  .0253 
0.0351 

0 .0395 
0  .0325 
0 .0419 

0.0456 
0.0384 
0  .0476 

0  .0510 
0.0435 


Standard  error  of  the  estimate  of  the  mean 
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TABLE  9 

Simulated  Means  and  Standard  Deviations  for  R(T; 
and  Their  Approximating  Values 


Nl  =  5'  N2  =  5'  Xl  =  10,  X2    =    40,  yl  =  25,  y2  =  125 


TIME  T 

0.01       Simulation 

CLT 
HT 
0.025      Simulation 

CLT 
HT 
0.0375     Simulation 

CLT 
HT 
0.0500     Simulation 

CLT 
HT 
0.0  625     Simulation 

CLT 
HT 


Mean 

0.0717 
(0.0001) * 

0  .0717 

0  .0720 

0.1788 
(0.0004) 

0  .1793 

0  .1799 

0.2690 
(0.0005) 

0 .2684 

0.2699 

0.3588 
(0.0006) 

0.3585 

0. 3599 

0.4481 
(0.0007) 

0.4481 

0  .4498 


Std.  Dev 
0.0134 

0.0237 
0.0133 
0.0280 

0.0375 
0 .0271 
0.0376 

0.0459 
0.0359 
0 .0456 

0.0530 
0.0433 
0.0527 

0 .0529 
0.0496 


Standard  error  of  the  estimate  of  the  mean 


TABLE  10 

Simulated  Means  and  Standard  Deviations  for  R(T1 
and  Their  Approximating  Values 


Nl  =  5'  N2  =  5'  Xl  =  40,  X2    =  10,  yl  =  125,  M2  =  25 


TIME  T 

0.01      Simulation 

CLT 
HT 
0.025     Simulation 

CLT 
HT 
0.0375    Simulation 

CLT 
HT 
0.0500    Simulation 

CLT 

HT 
0.0625    Simulation 

CLT 
HT 


Mean 

0.0725 
(0.0001) * 

0  .0724 

0.0720 

0.1816 
(0.0004) 

0.1810 

0.1799 

0 .2717 
(0.0005) 

0.2714 

0.2699 

0.3622 
(0.0006) 

0.3619 

0.3599 

0  .4522 
(0.0008) 

0.4524 

0.4498 


Std.  Dev. 
0.0138 

0.0249 
0  .0133 
0.0290 

0  .0394 
0.0271 
0  .0389 

0  .0482 
0  .0359 
0.0472 

0.0557 
0.0433 
0  .0557 

0.0623 
0  .0496 


Standard  error  of  the  estimate  of  the  mean 
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TABLE  11 

Simulated  Means  and  Standard  Deviations  for  R(T) 
and  Their  Approximating  Values 
for  the  One-Type  Model 


N  =  10,  A  =  25,  u  =  75 


TIME  T 


Mean 


Std.  Dev 


0.1 


0.025 


Simulation 

CLT 

HT 
Simulation 


CLT 
HT 
0.0375     Simulation 

CLT 

HT 
0.0500     Simulation 

CLT 

HT 
0.0625     Simulation 

CLT 

HT 


0.0701 
(0.0002) * 

0.0701 

0  .0700 

0.1766 
(0.0005) 

0.1752 

0.1750 

0  .2620 
(0.0008) 

0.2628 

0.2625 

0.3501 
(0.0009) 

0  .3504 

0  .3500 

0.4388 
(0.0011) 

0.4380 

0.4375 


0.0139 

0.0206 
0.0135 
0  .0263 

0.0325 
0.0258 
0.0361 

0.0398 
0.0330 
0.0430 

0.0460 
0.0390 
0.0478 

0.0514 
0.0441 


Standard  error  of  the  estimate  of  the  mean 


36 


To  assess  the  quality  of  the  normal  approximation  to  the 
distribution  of  R(T)  for  the  two-type  model,  the  a-quantiles 
for  each  two-type  approximating  normal  distribution  were 
computed.   The  relative  frequency  of  being  less  than  or  equal 
to  each  a-quantile  was  then  computed,  using  the  simulated 
data.   The  results  appear  in  Tables  12-15.   From  the  heavy 
traffic  approximation  to  the  mean  it  follows  that  approximately 
7  jobs  are  being  processed  with  the  tagged  job.   Thus,  all 
the  cases  considered  are  really  moderate  traffic  cases.   The 
tables  indicated  that  the  HT  approximation  tends  to  describe 
the  quantiles  better  than  does  the  CLT.   However,  as  is 
expected,  the  CLT  improves  with  larger  T. 
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TABLE    12 

Simulated   Probability    (Relative   Frequency)     that    the 
Response   Time    is    Less    than   or    Equal    to    the   Approximating 

ct-Quantiles 


Nl    =    5'    N2    =    5'    Xl    =    20'    A2    =    30,    yl    =    50,    y2    =    10° 
TIME    T  a:  0.10         0.25         0.50         0.75         0.90         0.95         0.99 


0.010                 CLT  0.041  0.137  0.433  0.838  0.998  1.0  1.0 

HT  0.109  0.221  0.443  0.719  0.915  0.977  1.0 

0.0250              CLT  0.062  0.170  0.425  0.765  0.969  0.997  1.0 

HT  0.103  0.223  0.438  0.708  0.911  0.975  1.0 

0.0375              CLT  0.079  0.191  0.429  0.754  0.952  0.991  1.0 

HT  0.118  0.229  0.447  0.713  0.905  0.969  0.999 

0.0500              CLT  0.075  0.181  0.417  0.724  0.939  0.986  1.0 

HT  0.106  0.222  0.431  0.693  0.897  0.965  0.998 

0.0625              CLT  0.081  0.198  0.436  0.745  0.932  0.980  1.0 

HT  0.115  0.238  0.451  0.719  0.898  0.960  0.997 


TABLE  13 

Simulated  Probability  (Relative  Frequency)  that  the 
Response  Time  is  Less  than  or  Equal  to  the  Approximating 

a-Quantiles 


Nl  =  5'  N2  =  5'  Al  =  30,  A2  =  20,  yl  =  100'  y2  =  50 


TIME  T    a:  0.10  0.25  0.50  0.75  0.90  0.95  0.99 

0.0100    CLT  0.049  0.159  0.457  0.869  0.998  1.0  1.0 

HT  0.114  0.229  0.442  0.718  0.919  0.979  1.0 

0.0250    CLT  0.078  0.194  0.462  0.809  0.978  0.997  1.0 

HT  0.113  0.228  0.442  0.717  0.916  0.973  0.999 

0.0375    CLT  0.088  0.204  0.455  0.787  0.962  0.995  1.0 

HT  0.111  0.224  0.432  0.702  0.909  0.965  0.999 

0.0500    CLT  0.088  0.212  0.459  0.771  0.955  0.989  1.0 

HT  0.111  0.226  0.434  0.696  0.898  0.961  0.998 

0.0625    CLT  0.093  0.212  0.465  0.773  0.945  0.986  1.0 

HT  0.110  0.222  0.432  0.705  0.889  0.954  0.997 
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TABLE  14 

Simulated  Probability  (Relative  Frequency)  that  the 
Response  Time  is  Less  than  or  Equal  to  the  Approximating 

a-Quantiles 


Nl  =  5'  N2  =  5'  Al  =  10'  X2  =  4°'  Ul  =  25/  U2  =  125 


TIME  T    a  0.10  0.25  0.50  0.75  0.90  0.95  0.99 

0.010  CLT  0.027  0.119  0.447  0.906  1.0  1.0  1.0 

HT  0.110  0.230  0.458  0.737  0.931  0.982  1.0 

0.0250    CLT  0.061  0.178  0.451  0.824  0.991  0.999  1.0 

HT  0.121  0.236  0.464  0.734  0.930  0.985  1.0 

0.0375         CLT  0.068  0.191  0.441  0.790  0.979  0.999  1.0 

HT  0.117  0.239  0.451  0.725  0.927  0.983  1.0 

0.0500         CLT  0.081  0.189  0.430  0.776  0.975  0.997  1.0 

HT  0.117  0.232  0.444  0.729  0.935  0.985  1.0 

0.0625         CLT  0.085  0.199  0.442  0.772  0.960  0.995  1.0 

HT  0.121  0.237  0.456  0.737  0.923  0.979  1.0 
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TABLE  15 

Simulated  Probability  (Relative  Frequency)  that  the 
Response  Time  is  Less  than  or  Equal  to  the  Approximating 

a-Quantiles 


Nl  =  5'  N2  =  5'  Al  =  40,  A2  =  10'  yl  =  125/  U2  =  25 


TIME  T    a:  0.10  0.25  0.50  0.75  0.90  0.95  0.99 

0.0100         CLT  0.024  0.116  0.445  0.907  1.0  1.0  1.0 

HT  0.110  0.228  0.433  0.710  0.905  0.973  1.0 

0.0250         CLT  0.054  0.164  0.435  0.820  0.993  1.0  1.0 

HT  0.112  0.219  0.422  0.682  0.901  0.967  1.0 

0.0375         CLT  0.070  0.186  0.428  0.795  0.981  0.999  1.0 

HT  0.113  0.224  0.414  0.693  0.899  0.967  0.999 

0.050  CLT  0.075  0.182  0.439  0.784  0.972  0.998  1.0 

HT  0.110  0.216  0.422  0.689  0.896  0.969  1.0 

0.0625         CLT  0.083  0.200  0.443  0.765  0.959  0.997  1.0 

HT  0.118  0.222  0.427  0.683  0.891  0.958  1.0 
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APPENDIX 

HEAVY  TRAFFIC  APPROXIMATION  BY  CONVERGENCE- 
OF- SEMI GROUPS  METHODOLOGY 


The  purpose  of  this  appendix  is  to  outline  a  mathematical 
framework   upon  which  the  heavy  traffic  approximations  of  this 
paper  may  be  rigorously  based.   The  approach  is  to  use  an 
analytical  theory  of  convergence  of  semigroups  of  operators 
apparently  first  applied  to  queueing  problems  by  Barman  (19  79) 
in  a  regrettably  unpublished  thesis.   See  also  Lehoczky  and 
Gaver  (1981)  where  the  technique  is  used  to  obtain  results 
concerning  a  data-voice  traffic  sharing  multichannel  system. 
The  theory  of  semigroups  of  operators  is  introduced  in  Feller 
(1971),  and  detailed  in  Dynkin  (1965);  the  convergence  ideas 
are  discussed  in  Trotter  (1974)  and  Kato  (.19  76)  .   The  basic 
notion  is  that  the  state  variable  of  a  process,  say  the  work 
time  process  of  Section  4.1,  iX.(w;N),  w    0),  is  one  of  a 
sequence  of  birth  and  death  Markov  processes  indexed  by  system 
size  N.   Given  such  a  sequence  of  Markov  processes,  <  '  X.  (w;N)  }>  , 
each  with  its  appropriate  state  space,  S  ,  it  is  desired  to 
show  that  the  corresponding  sequence  of  probability  transition 
functions  converges  to  that  of  some  limiting  process  that  has 
state  SDace  S  ;  in  the  present  case  S   =  TR  .   The  limiting 

OO  ~  CO  ■+■  3 

process  under  the  normalization  of  {X(w;N)}  chosen  will  be  a 
particular  diffusion  process,  namely,  in  the  present  heavy 
traffic  situation  the  multivariate  Ornstein-Uhlenbeck . 

The  Trotter-Kato  theory  of  convergence  deals  with  the 
convergence  of  infinitesimal  operators  A^  of  the  normalized 
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processes  {X^(w);N>:   if  a  f  -*■  A^f  in  sup  norm  for  a  suitable 

class  of  test  functions  (e.g.  f(z):  TR      -*■  TR   m-times  continuously 

dif  ferentiable ,  m  >_  3 ,  that  vanish  identically  outside  a  bounded 

subset  of  TR   and  further  such  that  the  functions,  f,  together  with  their 

first  and  second  derivatives  do  not  increase  faster  than  some 

fixed  power  of  z)  it  can  be  concluded  that  the  semigroups 

converge,  and  hence  the  Markov  probability  transition  functions 

themselves  converge. 

We  now  proceed  with  the  formal  calculation  of  the  limiting 
generator  for  our  normalized  process.   Invoke  (4.4)  so 

X. (w)  -  N.m.(w)      XN(w)  -  N£.m.(w) 

ZN(w)   =   -* i-i =   -1 *-* .  (a-1) 

1  /FT  /T~  /N 

l  l 

By  definition,  for  z  =  ( z, , z_ , . . . , z)  and  f  in  the  above  class 

1       2.  K 

Af(z)       =       lim{E[f (ZN(w+A) ) |ZN(w)    = z ]     -    f(z)}|    ■  (A-2) 

A-^0 

Given    Z. (w)    =    z. ,    and   C. (w,w+A)    represents    the   change    in    X. (w) , 

XN(w)     +    C. (w,w+A)     -    N.m. (w+A) 
N ,     ,  .  s  l  i  li 


Z     (w+A) 


l 


C     (w) 

— +    z.     -    /N~  m!  (w)A  +    o(A)     .  (A-3 

i 


Consequently,    for    z    such    that    z.     >    -/N~~ m.  (w)  A    - 
^  J  —  l  li 


L 


'N. 


i 
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K 


A^f (z)       =      lim[    I    {f  (z    ,. 
iN  A+0    i=l  x 


.,z.    + /FT/Nm!  (w)  A,.  .  .  ,zTJ  A.  (N)  A 

1         /I7/N  x        x  K      x 


+    f (z     , .  .  .  ,  z  . 

1  X       /T~/N 


/£./Nm!  (w)  A,  .  .  .  ,  z__)  y  .  (N)  A} 
11  K       1 


+    f  (z    -/17/Nm'  (w)  A,  .  .  .  ,  z    -/TT/Nni'  (w)  A) 
111  K  K  K 


K 


1    -      £    {(A. (N)    +   y. (N) )A} 
i=l        x  1 


-    f  (z1,z2,  .  .  .  ,  z    )    +    o(A)  ]-£ 


:a-4 


where    for    simplicity    (and   generality)    we    abbreviate 


Ai(N! 


K 

=      A  .  [N.  -  n.  m.  (w)-/N.z.J  (    J     (N  .m  .  (w) +/N  .  z  .  ) 
iiii  l  l      .^       ]    j  J    J 


(A-5) 


(A .N)  [I.  (1-m.  (w)  ) 


/I  .  K 

— -  z. ]  (N      7    U .m.  (w)    + 
/N       x  j=l      3D 


J. 


z  . 


/N        ^ 


A-6) 


and 


y . (N)       =       y . (N.m. (w)     +    /Nz . 
l  ill  l 


yiN(£imi(w)    + 


l 


/N        X 


Z.)      .  (A-7) 


Upon  passage  to  the  limit  via  Taylor  series  expansion  it  is 
seen  that 
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K 
A^f(z)   =    I    {f(z ,...z   +  — ± ,...,z    )A. 


n: 


+  f (z. , . . . ,z.  - 

1      1 


,  .  .  .  , z  ) y  .  (N)  } 

/£~/N       K   X 


+  f (z   ,  .  .  .Z  .  ,  . 


,zk,[- 


K 

I    (X  .  (N)  +  y .  (N) )  ] 
=  1   3        D 


K 
•-1    Z  • 


£   /N  m' (w' 
J      J 


A- 


Note  that  no  specific  normalization  has  been  utilized  up  to 
this  point.   Now,  however,  invoke  the  HTN  of  (4.6)  and  utilize 
(A-6)  and  (A-7)  specifically;  allowing  N  to  become  large, 


*N 


f  (z 


K       A  .  (N)  -  y  .  (N)     ,         A  .  (N)  +  y  .  (N 

1=1    l   /£ .  /N  l 

l 


-  f '   /L  /Nm!  (w)  } 
z  .     1       1 


(A-9 


K 


K 


y   f  {/N[A.V£.  (1-m.  (w)  )   j      £.m.(w)  -  y  .  /FTm.  (w) 
i=l   zi      x   x     x      j  =  l   3  3        X   X  x 

K  K 

-/Z.m!(w)]  +  A  !/T~(l-m.  (w))(  7  /£.z.)-A!z.      Jl.m.(w) 

-  y.z.  +  0(N_1/2)}  +  i   I   f"  (AMl-m.(w)) 
li  Z.-.Z.1     l 


K 


-1/2 
x   y  9.    m  (w)  +  y  .m.  (w)  +  0(N   '  )  } 

D=i    3  j  1  1 


v. 


Choose  the  functions  m.  so  that  they  satisfy  the  system  of 
differential  equations  (4.5)  .   Let  N  ■*  °° .   Then  for  f  in  the 
above  class  of  functions,  the  operator  A^  converges  to  yield 

K  K  K 

A  f(z)   =      f  {A.'/JT  (1-m.  (w)  )      /T7z.  -  A!z.      I  .m.  (w)  -y  .  z .  J 

i=l   zi   x   1     x     i=l    J    3  x  x  i=i   3  3      ii 


K  K 

+  ^  f"  (A!  (1-m.  (w))  (  y  Jl.m.(w))  +  y.m.(w)} 

2   ^ ,  z  .   l     l        i   1  J         ii 

1  =  1  l              j=l   J  J 


The  operator  A^  is  the  infinitesimal  operator  of  the  diffusion 
whose  stochastic  differential  equation  is  (4.6)  [cf.  Arnold 
(1974),  page  152].   The  Trotter-Kato  theorem  can  now  be  applied 
to  assert  that  the  semigroups  converge  [cf.  Burman  (1979)]  . 
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